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ABSTRACT 

We have constructed a detailed radiative transfer disk model which repro- 
duces the main features of the spectrum of the outbursting young stellar object 
FU Orionis from ~ 4000 A to ~ 8jj,m. Using an estimated visual extinction 
Ay ~ 1.5, a steady disk model with a central star mass ~ O.SMq and a mass ac- 
cretion rate ~ 2 X 10-^MQyr-\ we can reproduce the spectral energy distribution 
of FU Ori quite well. Higher values of extinction used in previous analysis (Av/~ 
2.1) result in spectral energy distributions which are less well- fitted by a steady 
disk model, but might be explained by extra energy dissipation of the boundary 
layer in the inner disk. With the mid-infrared spectrum obtained by the Infrared 
Spectrograph (IRS) on board the Spitzer Space Telescope, we estimate that the 
outer radius of the hot, rapidly accreting inner disk is ~ 1 AU using disk models 
truncated at this outer radius. Inclusion of radiation from a cooler irradiated 
outer disk might reduce the outer limit of the hot inner disk to ~ 0.5 AU. In 
either case, the radius is inconsistent with a pure thermal instability model for 
the outburst. Our radiative transfer model implies that the central disk temper- 
ature Tc > 1000 K out to ~ 0.5 — 1 AU, suggesting that the magnetorotational 
instability can be supported out to that distance. Assuming that the ~ 100 yr 
decay timescale in brightness of FU Ori represents the viscous timescale of the 
hot inner disk, we estimate the viscosity parameter to be a ~ 0.2 — 0.02 in the 
outburst state, consistent with numerical simulations of the magnetorotational 
instability in disks. The radial extent of the high M region is inconsistent with 
the model of Bell & Lin, but may be consistent with theories incorporating both 
gravitational and magnetorotational instabilities. 
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Introduction 



The FU Orionis systems are a small but remarkable class o f variable you np; objects 
which undergo outbursts in optical light of 5 magnitudes or more (iHerbigl 119771 ). with a F- 
G supergiant optical spectra and K-M supergiant near-infrared spectra dominated by deep 
CO overtone absorption. While the rise times for outbursts are usually very short (~ 1-10 
yr), the decay timescales range from decades to centuries. The FU Ori objects also show 
distinctive reflection nebulae, large infrared excesses of radi ation, wavelength dependen t 
spectral types, and "double-peaked" absorption line profiles ( iHartmann fc Kenyoru Il996l ). 
The frequency of these outbursts is uncertain; in recent years an increasing number of heavily 
extincted potential FU Ori objects have been identified on the b a sis of their spectroscopic 
characteristics at nea. r -infrared wavelengths (IKenyon et al.l Il993l : iReipurth fc AspinI 119971 : 
Sandell fc A^ll998l : kspin fc Reipurthlboosh . 



The accretion disk model for FU O ri objects proposed by Hartmann & Kenyon (1985, 
1987a, 1987b) and iKenyon et al.l (119881 ) can explain the peculiarities enumerated above in 
a straightforward manner. Outbursts are known in other accreting disk systems and may 
be the result of a common mechanism (e.g.. Bell & Lin 1994). The high temperature of the 
inner disk produces the observed F-G supergiant optical spectrum, while the cooler outer 
disk produces an infrared spectrum having the spectral type of a K-M supergiant. The 
Keplerian rotation of the disk can produce double-peaked line profiles as often observed, 
with peak separation decreasing with increasing wavelength of observations, since the inner 
hotter disk which produces the optical spectrum rotates faster th an the outer cooler disk 
which produces the infrared spectrum (IHartmann fc Kenyonlll996l ). 



While the accretion disk model has been successful so far, it is important to continue 
to test it and to derive further insights into the origins of the accretion outbursts as new 
observations become available. In particular, the mid-infrared spectrum obtained with the 
Infrared Spectrograph (IRS) on the Spit zer Space Telescop e provides important constraints 
on the outer edge of the hot inner disk (j Green et al.l |2006| ) , which in turn can test theories 
of outbursts. For example, the thermal instability model of iBell fc LinI (119941 ) predicts an 
ou ter radius for the o utburst region of only ~ 2OR0, which is small er than that est i mated 
by iGreen et al.l (120061 ) on the basis of simple blackbody disk models. lArmitage et al.l (1200 ll ) 
suggested that a combination of gravitational instability (GI) with triggering of the mag- 
netorotational instability (MRI) might also explain FU Ori outbursts; in this model the 
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high accretion rate region in outburst would be much larger, of order 0.5 AU. Vorobyov & 
Basu (2005,2006) suggested that FU Ori outbursts could also be produced by the accretion 
of clumps formed in a gravitationally unstable disk. This model predicts high accretion 
rates at its inner radius of 10 AU. These differing predictions for the size of the outburst 
region should be testable with detailed models of the spectral energy distribution (SED). 
Furthermore, if we can constrain the size of the outburst region, the timescales of decay 
may provide quantitative estimates of the viscosity, vital for understanding the evolution of 
young protoplanetary disks. 

In this paper we present new steady, optically-thick accretion disk spectrum calculations, 
and compare the results with observations from the ultraviolet to the mid-infrared of FU 
Ori— the prototype of these systems. In §2 we describe the observational material used to 
constrain the models, while in §3 we describe the methods used to calculate the disk models. 
We show in §4 that the steady disk model reproduces the variation of the observed spectral 
features over this very large wavelength range quite well. In §5 we consider some implications 
of our results for disk viscosities, outburst mechanisms, and disk masses, and summarize our 
conclusions in 56. 



2. Optical and infrared data 



An essential part of the spectral energy distribution (SED) of FU Ori for our model is 
the IRS spectrum from Spitzer observed on March 4, 2004 (see Green et al. 2006 for details). 
As FU Ori has been fading slowly over the last 60 years, we need to obtain data at other 
wavelengths near the same time as the IRS spectrum to construct a complete spectral energy 
distribution (SED). 



Optical photometry (UBVR) was obtained in 2004 at the Maidanak Observatory (jibrahimov 



19991 : iGreen et al.l 120061 ). As FU Ori exhibits small, irregular variability on timescales much 
less than a year, we averaged the Maidanak data to form the mean values (with RMS un- 
certainty in each band ~ 0.02 magnitude) used in this paper. In addition, we include an 
optical spectrum of FU Ori obtained on March 27th, 2004, near the time of the IRS obser- 
vations, using the 1.5m telescope of the Whipple Observatory with the FAST Spectrograph 
(IFabricant et al.lll998l ). equipped with the Loral 512 x 2688 CCD. This instrument provides 
3400 A of spectral coverage centered at 5500 A, with a resolution of 6 A. The spectrum 
was wavelength calibrated and combined using standard IRAF routines 0. The spectrum 



^IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the As- 
sociation of Universities for Research in Astronomy, Inc., under cooperative agreement with the National 
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was corrected for the relative system response using the IRAF sensfunc task. The optical 
spectrum was placed on an absolute flux scale using the averaged 2004 UBVR photometry. 

Near-infrared fluxes were determined from several sources. 2MASS JHK photometry is 
avail able from October 7th, 1 999, five years before IRS observations. We also include data 
from lReipurth fc AspinI (120041 ) on December 15, 2003, using the NAO J 8 m Subaru Telescope 
on Mauna Kea and the Infrared Camera and Spectrograph (IRCS). Finally, we obtained JHK 
magnitudes on December 17th 2005, using the Near-IR Imager/Spectrometer TIFKAM on 
the 2.4 m MDM telescope at Kitt Peak. Although the 2.4m telescope had to be defocused 
to avoid saturation because FU Ori is so bright, our JHK magnitudes of 6.57, 5.91 and 5.33, 
respectively, are in reasonable agreement with the JHK^ 2MASS magnitudes of 6.519, 5.699, 
5.159, respectively, indica ting that the FU Ori flux es have not changed dramatically over 
that six year period. The iReipurth fc AspinI (120041 ) observations are on a different system, 
but when converted to fluxes yield reasonably similar results. We therefore adopted the 
2MASS magnitudes as our standard in this wavelength region. 

We then incorpo rated the KSPEC spect rum (resolution R~500, 1.15-2.42/im ) observed 
on Dec 15th, 1994 by I Greene fc Ladal (119961 ) and the observation with SpeX (2.1-4.8/im) on 
the Infrared Telescope Facility (IRTF; Rayner et al. 1998) obtained during the nights of 
January 5-6th, 2001 (see MuzeroUe et al. 2003 for details). We placed these spectra on an 
absolute scale by convo lving the spectra w ith the relative spectral response curve (RSRs) 
of the 2MASS system JCohen et all Wm and tying the overall resuh to the 2MASS K^ 
magnitude. 

The collected observations are shown in Figured! where we assume an Ay of 1.5 instead 
of the usually assumed value of ~ 2.2 (Kenyon et al. 1988; KHH88). The previously discussed 
variation of spectral type with wavelength is immediately apparent in Figure (H the optical 
spectrum corresponds to a G type spectrum as estimated from the atomic metal lines while 
the near-IR and IRS wavelength ranges are dominated by molecular bands indicative of a 
much later spectral type. While the optical and infrared spectra are dominated by absorption 
features, dust emission dominates longward of 8 /xm, as shown by the silicate emission 
features at 10 and 20 fim. 
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Model 



In this paper we only model the absorption spectrum (A<8 /xm) which is produced by the 
centrally-heated accretion disk. The silicate emission features are produced by externally- 
heated regions of the outer disk or possibly a circumstellar dusty envelope (Kenyon & Hart- 
mann 1991; Green et al. 2006); we defer modeling of this region to a subsequent paper. 

We follow the method of Calvet et al. (1991a,b) to calculate the disk spectrum. In 
summary, we calculate the emission from the atmosphere of a viscous, geometrically thin, 
optically thick accretion disk with constant mass accretion rate M around a star with mass 
M and radius R. The disk height H is assumed to vary with the distance from the rotational 
axis of the star R as H = Hq{R/ RiY^^ , where we take Hq = O.li?^ and Ri is the radius of 
the central star. This approximation is not very accurate but it only affects the local surface 
gravity of the disk atmosphere, which has only a small effect on the emergent spectrum. 
Viscous dis sipation in the atmosphere is neglected, which is a good approximation for cases 
of interest (IHartmann fc Kenyonlll99ll ). Thus, radiative equilibrium holds in the disk atmo- 
sphere, and the surface flux is determined by the viscous energy generation in the deeper 
disk layers. This constant radiative flux through the disk atmosphere can be characterized 
by the effective temperature distribution of the steady optically-thick disk. 



. _ 3GMM 
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This equation predicts that the maximum disk temperature Tmax occurs at 1.36-Ri and then 
decreases to zero at R = R^. The latter is unphysical, and so we modify equation ([1]) so 
that when the radius is smaller than 1.36 Rj, we assume that the temperature is constant 
and equal to T = T^ax- (e.g., KIIH88) The vertical temperature structure at each radius is 
calculated using the gray- atmosphere approximation in the Eddington limit, adopting the 
Rosseland mean optical depth r. 

The opacity of atomic and mol ecular lines has been cal c ulated using the Opacity Dis- 
tribu tion Function (ODF) method (jCastelh fc Kuruczl 12004 : Isbordone et all 12004 : ICasteUi 
20051 ). Briefly, the ODF method is a statistical approach to handling line bla nketing when 
millions of lines are present in a short wavelength range (IKurucz et al.lll974l ). For a given 
temperature and pressure, the absorption coefficient for each line is exactly computed, then 
the profiles of all the lines in each small interval Au are rearranged monotonically as a func- 
tion of u. The opacity increases as the frequency increases in this Au. A step function with 
12 subintervals in frequency is used to represent this monotonic function. The height of 
every step is the averaged opacity of the monotonic function in this subinterval. The width 
of the 12 subintervals are Au/lO for 9 intervals and Az//20, Au/30, Au/QO for the last 3 
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intervals where the monotonic function increases steeply and reaches the maximum. Thus, 
in every interval Az/ we obtain 12 representative averaged opacities at a given temperature 
and pressure. We then construct an opacity table as a function of A for each tem perature 
and p ressure. The wavelength grid is the same as used in the code DFSYNTHE ( ICastelli 
2005h from 8.976 nm to 10000 nm with 328 BIG intervals (resolution 30-100) and 1212 LIT- 
TLE intervals (resolution 60-500). However, we extend the temperature grid and pressure 
grid to lower temperatures and pressures (logT ranging from 1.5 to 5.3, logP from -24 to 
8 in cgs units) than DFSYNTHE, so that it can be used in our disk models for which the 
temperatures and pressures reach lower values than in typical stellar atmosphere models. 
Because the line opacity usually varies by several orders of magnitude within a line width, 
the ODF method is substantially more accurate than either straight or harmonic means. 



The line hst is taken from Kurucz's CD-ROMs Nos. 1, 15, 24, and 26 flKurucd 120051 ). 
Not only atomic lines but also many molecular lines (C2, CN, CO, H2, CH, NH, OH, MgH, 
SiH, SiO, TiO, H2O) are included . The opacities of Ti O and H2O , the mos t impo rtant 
molecules in the infrared, are from [Partridge fc Schwenkd (119971 ) and ISchwenkd (119981 ). We 
do not include a detailed calculation of the dust condensation process; instead, we add the 
dust opacity when T<1500 K (the cond ensation temp e rature of si licates at typical disk 
densities). We use the ISM dust model of iDraine fc Led (Il984j . 119871 ) to represent the dust 
in the disk. This assumption is unlikely to be correct but mostly affects the continuum at 
the very longest wavelengths and in particular the silicate emission features, which we do 
not model in this paper. 

At low temperatures complex chemical processes occur which are not included in the 
Kurucz data. We have not calculated the low temperature molecular opacity in detail; 
instead, we have assumed that the abundance ratio between different types of molecules 
below 700 K is the same as the ratio at 700 K. This is unimportant for our purposes because 
dust opacity dominates at such low temperatures. 

In Figure [2] we show the Rosseland mean opacity calculated for solar abundances 
and a turbulent veloci t y v = 2kms~^ for different values of p/T| to be compared with 
Alexander fc Fergusonl (jl994j ). where p is the density in gcm~^ and Tg is the tempera- 
ture in millions of deg r ees. The results closely match t he more detailed calculations by 
Alexander fc Fergusonl (119941 ) and iFerguson et al.l (120051 ) except near 1500 K, where our 
dust opacity rises more rapidly due to our neglect of dust condensation processes. 



Calvet et al.l (jl991bl ) showed that the near-infrared spectrum of FU Ori could be well- 



modeled with the then-current water vapor opacities and a similar calculational method to 
the one we are using here. We now have a much better set of opacities, and can treat the 
optical spectrum. To test the code at optical wavelengths, we compare the model spectrum 
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from the an nulus with T^f /■=53 00 K with the observed spectrum of SAO 21446 which is a Gl 
supergiant (jjacoby et aLlll984l ). The spectrum of SAO 21446 has been obtai ned using the 



Intens ified Reticon Scanner(IRS) on the No.l 90 cm telescope at Kitt Peak by lJacoby et al. 



( Il984l ). Here we have convolved this spectrum to a resolution of 300 to agree with our model 
spectrum. As shown in Figure [3l our model reproduces the continuum spectrum fairly well 
except for a few strong features. These discrepancies are probably due to the limitation of 
the grey atmospheric assumption. 

Figure H] shows the emergent intensities with 55° emergent angle from different annuli of 
our final disk model (parameters can be found in §4), which has 45 annuli . The 55° inclination 
angle is estimated from the near-infrared interferometric observations by lMalbet et al.l (120051 ) 
and mid-infrared interferometric observations by Quanz et al. (2006). The radii of these 
annuli are chosen to increase exponentially from r = Ri to 1000-Rj. At the first annulus 
(r = IRi), we can see a significant Balmer jump at 0.36 fim, similar to the spectra of early 
type stars. As the radius increases, the temperatures of the annuli become lower and lower. 
Infrared molecular features appear at larger annuli with effective temperatures less than 
5000K (r > 3Ri). For the outermost annuli with effective temperature less than 1500 K 
(r > 15i?j), the molecular features have almost disappeared because dust opacity dominates. 
The final spectrum is the result of the addition of the fluxes from each of these annuli 
weighted by the appropriate annular surface area. 

In addition to these calculations aimed at comparing with low-resolution spectra, we also 
calculated high-resolution spectra in restricted wavelength regions. For these calculations, 
the vertical structure of every annulu s is calculated as above, and then th is structure is 
imported into the program SYNTHE JXurucz fc Avrettlll98ll : lKuruczlll993[ l. SYNTHE is 
a suite of programs which solve the radiative transfer equations in LTE with a very high 
spectrum resolution (~500,000). Each annular spectrum is broadened with the rotational 
profile appropriate to a Keplerian disk and then the summed spectrum is calculated (e.g., 
KHH88). These calculations are very time-consuming, and were used only to produce a high 
resolution spectrum of a small wavelength range around 6170 A for the purpose of estimating 
the central mass, and for examining the CO first-overtone lines. 



4. Comparison with observations 

The SED of a steady, optically-thick disk model is determined by two parameters: the 
product of the mass accretion rate and the mass of the central star, MM, and the inner 
radius Ri. One observational constraint comes from the observed spectral lines and/or the 
peak of the SED (after reddening correction), which determines the maximum temperature 
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of the steady disk model, 

/ 3GMM \ 

The other observational constraint is given by the true luminosity of a flat disk L^, 

. F GMM 
cos z 2Ri 

where d is the distance of FU Ori, i is the inclination angle of the disk to the line of sight, 
and F is the observed total flux corrected for extinction. If the distance and inclination are 
known, we can solve for the inner disk radius and thus for the product MM. If we then 
further use the observed rotational velocities, it is then possible to derive independent values 
of the central mass M and the accretion rate M, in the manner outlined by KHH88. 

The best constraints on Tmax come from the optical spectrum, as the uncertain extinction 
makes it difficult to constrain by overall SED fitting. From the lines of our optical FAST 
spectr um we derive a spectral type of ~ G2 for FU Ori us ing the methods of iHernandez et al. 



( 120041 ). which agrees with previous results (jHerbig||l977l ). However, it is not straightforward 



to derive a value of extinction and thus Tmax from this spectral type determination; the disk 
spectrum is not that of a single (standard) star with a well-defined effective temperature, but 
instead it is a combination of hotter and cooler regions, with an overall spectrum that clearly 
varies with the wavelength of observation in a complicated way. To address this problem, 
we constructed several disk models with a modest range of Tmax and compared them with 
the observed optical spectrum dereddened by differing amounts of visual extinction Ay, as 
shown in Figure O 

The uppermost model in Figure [5] has Tmax = 7240 K, comparable to that of the FU 
Ori disk model of KHH88. The model optical spectrum is too flat longward of ~ 3900 A 
in comparison with the observations dereddened by the Ay = 2.2, suggested by KHH88, or 
for Ay = 1.9. At the other extreme, we flnd that the bottom model, with Tmax = 5840 K, 
matches the observed spectrum dereddened by only Ay = 1.3 quite well. However, we 
suspect that our gray-atmosphere approximation somewhat underestimates the amount of 
line blanketing in the blue optical region, especially shortward of the Ca II resonance lines 
at 3933 and 3968 A leading us to suspect that this agreement is somewhat misleading. We 
therefore provide our best estimate of the extinction as Ay = 1.5 ±0.2, and adopt the model 
which best matches the mean extinction with Tmax = 6420 K as our standard model. 

In this connection it is worth noting that KHH88 did not actually flnd a good flt 
for the optical spectrum of FU Ori using a GOI standard star and extinctions Ay > 2. 
They speculated that the discrepancy might be alleviated by an improved treatment of limb 
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darkening in the disk. However, Hartmann & Kenyon (1985) had found that a G2I standard 
star and Ay = 1.55 provided a good match to the optical spectrum between 3900 A and 
7400 A , consistent with our resuhs here. 

The accretion disk luminosity L depends upon the distance, inclination, and dereddened 
flux. We adopt a distance of ~ 500 pc, consistent with membership in the general Orion 
region (Herbig 1977). The inclination 55° is adopted (Malbet et al.2005, Quanz et al. 2006). 
Using these parameters, observed total flux and Ay = 1.5, we obtain the true luminosity 
Ld = 8.66 X lO^^ergss~^~226L0 according to equation [3l As mentioned before, fitting the 
observed spectrum yields T^ax = 6420 K. We derive MM = 7.2 x 10-^M|yr-i and Ri = 5Rq 
by solving equation [2] and equation [3] simultaneously. 

Figure [6] shows the predicted spectrum for the adopted parameters and several values of 
the outer radius. Our computed spectrum matches the observations at wavelengths short- 
ward of about 4/im, indicating that the innermost region of the inner disk is reasonably well 
reproduced by the adopted MM and Ri parameters. 

Because the optical and near-infrared SED is well-matched by the model, we may then 
proceed to estimate the central stellar mass and thus the accretion rate. We used the 
observed optical rotation from KHH88 to constrain the central star mass. We computed 
the line profiles around 6170 A using SYNTHE, cross correlated it with a non-rotating disk 
spectrum, and compared it with the observed cross-correlation line profile given in KHH88. 
The uncertainty in fitting the line profile widths is about ± 20%. However, systematic errors 
in the inclination angle are probably more important for the mass estimate, and for that 
reason we give an uncertainty O.3M0±O.1M0 ( consistent with the earlier estimate of KHH88 
assuming i = 50°). For a mass of 0.3 Mq, M ~ 2.4 x lO~^M0yr~^. Parameters of the model 
are given in Table 1. 

As shown in Figure [H we find that regardless of extinction and disk parameters, steady 
accretion models which fit the optical to near-infrared region predict too much emission in 
the IRS range for very large outer radii. This can oi ily be remedied b y truncating the hot 



inner disk. The estimate of the outburst models of iBell &: Lid (119941 ) suggested an outer 
radius of the high state of ~ 2OR0 (with somewhat different parameters for the accretion rate 
and inner disk radius); as shown in Figure [6] this truncation radius fails to explain the SED. 
A truncation radius of Rout ~ 210Rq ~ 1 AU provides a better fit to the flux at ~ 5 — 8/xm, 
although is s omew hat low in the 3/im band. Our results are reasonably consistent with 



Green et al.l (120061 ). considering the limitations of their blackbody modeling only out to 5 



microns. 



As discussed above, models with outer radii larger than ~2OOi?0 overpredict the flux at 
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~ 5 — 8/im. At even longer wavelengths, absorption features are no longer present; instead, 
sili cate dust en i ission features are seen in IRS spectra at ~ 10 and 20yum. As discussed 



by iGreen et al.l (120061 ) . these silicate features are signatures of heating from above, rather 
than from internal viscous dissipation, and may be produced in upper atmosphere of the 
outer disk which absorb light from the central disk. An outer disk is expected to be present, 
since a disk radius of ~ 0.5 — 1 AU is extremely small by standards of typical T Tauri 
disks, for which Rout ~ 100 AU or more are common (e.g., Simon, Dutrey, & Guilloteau 
2000). Moreover, such a small radius would imply a low- mass optically-thick disk, making 
it difficult to explain the observed submm flux from FU Ori (Sandell & Weintraub 2001). 
So, it is very likely that an outer disk is present with a lower accretion rate, which could 
contribute some flux in the 5 — 8/im wavelength range and thus reduce the radius of the hot, 
high-accretion rate region. An outer disk accretion rate of about 1O~^M0 yr~^ is comparable 
to some of the highest infall rates estimated for (not heavily embedded) protostars (e.g., 
Kenyon, Calvet & Hartmann 1993); higher accretion rates would imply implausibly short 
evolutionary timescales M/M ~ 3 x 10^ yr. In any event, if the outer disk accretion rate 
is not significantly smaller than the inner disk accretion rate, the problem of a large outer 
radius yielding too much fiux remains. More likely, the emission from the outer disk is 
dominated by irrad iation from the inner disk independent of the (lower) accretion rate there 



(ITurner et al.lll997l ). as indicated by the presence of emission rather than absorption features. 



In the iBell &: LinI (119941 ) thermal instability models for FU Ori outbursts, the outer disk 
accretion rate is of order 10^^ of the inner disk accretion rate in outburst. We have therefore 
investigated the effect of an outer disk accretion rate of 2.4 x 10~^Mq yr~^ on the SED. The 
dotted curves in Figure [6] show the results of adding in such outer disk emission from disk 
internal viscous dissipation. The effect is relatively small so that our estimate of the outer 
radius of the hot disk with the high accretion rate is relatively robust. Of course we cannot 
rule out other kinds of temperature distributions, such as a smooth decline of accretion rate 
from, say, 100 to 200 and the irradiation effect from the inner disk to the outer disk. We 
are presently carrying out more detailed calculations of the outer disk, including irradiation 
by the inner disk. Preliminary results indicate that the addition of the irradiated outer disk 
emission will decrease our hot disk outer radius estimate by no more than a factor of two 
(Zhu et al., in preparation). 

We explored the effects of a larger visual extinction by scaling the observations to 
Ay=2.2. As shown in Figure[71 this makes only a very slight difference to the long-wavelength 
spectrum and does not change our estimate of the outer radius. 



To examine our consistency with previous results, we calculated the CO first overtone 
high resolution spectrum with Rout = 210Rq using SYNTHE program and compared it with 
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the high resolution spectrum from iHartmann et al.l (12004 ). The comparison is shown in 
Figure [HI We found that we needed to adopt a turbulent velocity of 4 kms~^ to obtain lines 
that are deep enough in comparison with the observations. This is the sound speed of gas 
with temperature around 4000 K, slightly supersonic for the annulus with radius larger than 
4 stellar radius; Hartmann, Hinkle, & Calvet (2004) similarly found that slightly supersonic 
turbulence was needed to explain the first-overtone CO lines, and pointed out that this would 
not be surprising in the context of turbulence driven by the magnetorotational instability. 

In Figure [9] we show an expanded view of the wavelength region between ~ 3 and 
8/im. The model accounts for the overall shape of the SED reasonably well, and yields 
a reasonable strength for the 6.8/xm water vapor absorption feature. Howeve r, the 5.8/xm 



water vapor feature of the model is not as strong as it is in the observations ( (iGreen et al. 



20061 )). In this wavelength region the contributions from the dust-dominated regions of 



the disk are significant, and the contribution of this dust (relatively featureless) continuum 
emission could reduce the strength of the absor ption feature. To investiga te this possibility, 
we reduced the dust opacity to only 1% that of iDraine fc Led (Il984j . 119871 ) and recalculated 
the model (with a slightly larger outer radius). As can be seen in Figure [9l this reduction 
in dust opacity strengthens the absorption features, although the predicted 5.8/im feature 
is still not strong enough, perhaps suggesting some difficulty with the opacities. It is not 
clear whether we require dust depletion, because we have not considered the temperature 
dependence of dust condensation in our model. We also cannot rule out uncertainties due 
to our simple treatment of the vertical temperature structure of the disk. 



5. Discussion 

5.1. The steady accretion disk model and the inner boundary condition 

We have shown that the steady accretion disk model can quantitatively account for the 
observed SED of FU Ori over a factor of nearly 20 in wavelength, from 4000 A to about 
8/xm. This is accomplished with only three adjustable parameters: Ri, MM, and Rout- No 
alternative model proposed for FU Ori has explained the SED over such a large wavelength 
range. The accretion disk explanation is bolstered by this work. 

Within the context of the steady disk model, however, there remains the question of the 
appropriate inner boundary condition. The steady disk model considered here only accounts 
for half of the accretion lumino sity (L = GMM / {2R j)) potentially available from accretion 



onto a slowly rotating star. As iKenyon et al.l (119891 ) showed from ultraviolet spectra, there 



is no evidence in FU Ori for boundary layer emission which would account for this missing 
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luminosity, up to ~ 10-20% of the total luminosity; and if our lower estimate of extinction is 
correct, the potential boundary layer emission would be much less. This leaves as possibilities 
that the missing accretion energy is going into spinning up the central star, expanding the 
central star, being radiated over a larger radial distance of the disk, or some combination of 
all three effects. 

Hartmann & Kenyon (1985; also Kenyon et al. 1989) pointed out that the 5 inner 
radius implied by steady disk models is considerably larger than that of the typical T Tauri 
star {R ~ 2R0), and suggested that this might be due to the disk dumping large amounts 
of thermal energy into the star, expanding its outer layers. 

Popham et al. (1993, 1996) suggested that the the missing boundary layer energy could 
be distributed over a significant range of radii in the disk at high accretion rates. They 
found that as M increases to ~ 10~^ MQyr"^ the dynamical boundary layer (the radial 
extent of the region where Q drops from Keplerian angular velocity to 0) grows to 10-20 % 
of a stellar radius and the thermal boundary layer (the radial width over which the boundary 
layer luminosity is radiated) grows to the point that it becomes impossible to distinguish 
the boundary layer from the inner parts of the disk. Popham et al. (1996) showed that 
this model would imply somewhat slower than Keplerian rotation in the innermost disk. In 
this paper we find no evidence for extra radiation if Av'=1.5. We cannot rule out a slightly 
larger Ay, which would allow some extra dissipation of kinetic energy in the inner disk. 
However, the shape of the spectrum in the blue-optical region (Figure ED suggests that the 
higher extinction value is not preferred, limiting the amount of excess radiation to a modest 
fraction of the total luminosity. We speculate that the missing energy is mostly going into 
expanding the outer stellar layers, with perhaps a small amount of excess radiation that is 
difficult to discern given the uncertainties in the extinction. 



5.2. Rr„,t and a 



Our estimate of the outer radius of the high state of FU Ori is around 200 R©, although 
the irradiation effe ct from the hot in ner disk to the outer disk may decrease this value to no 



less than 100 Kq. iBell fc Liru (119941 ) assume low viscosity parameters (a~10 ^) and derive 



a small outer radius ( Ro„4~2OR0) for the high M, hot inner disk during outbursts. Our 
model of FU Ori is inconsistent with such a small hot region; instead, we require the high 
M region to be an order of magnitude larger in radius. 



FU Ori has been fading slowly in brightness over the last ~ 70 years (jlbrahimovlll999l ). 
If we attribute this fading to the emptying out of mass from the inner hot disk onto the 
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central star, this timescale gives us a way to constrain the viscosity in this region. We make 
the simple assumption that the decay timescale is the viscous timescale 

~ RVu , (4) 

where the viscosity is z/ = acl/Q, a is the viscosity parameter (Shakura & Sunyaev 1973), 
Cs is the sound speed in the disk midplane, and Q is the Keplerian angular frequency at 
R. We use the isothermal sound speed with a mean molecular weight of 2.3, appropriate 
for molecular hydrogen. We then need an estimate of the central temperature to evaluate 
ty. Our radiative transfer model shows that the effective temperature is Tg// ~ 800 K at 
the outer radius of i? ~ 2OOR0. This is a lower limit for the central temperature as the 
absorption features in the spectrum show that the central temperature is higher than the 
surface temperature. The disk is likely to be very optically thick in this region, such that the 
central temperature is considerably larger than the effective temperature due to radiative 
trapping. We choose 1500 K as the central temperature above which MRI is active and thus 
high a can sustain high mass accretion rate. We therefore find 

^ '''''' " (o^) ' (21^) ' (15^) (iFt)" ■ 

where M is the mass of the central star, R is the outer radius of the high mass accretion 
disk, T is the central disk temperature of the high mass accretion disk, and a is the viscosity 
parameter. If we estimate a typical decay timescale of FU Ori as ~ 100 yr, then we find 
a ~ 0.14. For the smaller hot disk radius with an irradiated outer disk, our preliminary 
results yield (Zhu et al., in preparation) R >0.5 AU and thus a ~ 0.1. 

It is important to recognize that equation ([5]) is only an order of magnitude estimate. 
The decay timescale of accretion can be affected by whatever mechanism makes the transition 
from high to low state; for example, in the thermal instability model, the decay time may 
not be the timescale for emptying out half the material in the affected region of the disk, 
but only that (smaller) amount necessary to shut off the thermal instability. N evertheless, 



i t is su ggestive that we find a larger value of a in the high state than the 10 of iBell &: Lin 



( 119941 ). 



A value of a ~ 0.2 — 0.02 is roug hly consistent with estimates from compact systems 



with accretion disks (IKing et al.l 120071 ) as well as with numerical simulations of the magne- 
torotational instability (MRI) in ionized disks (Balbus & Hawley 1988). In this connection, 
we note that the requirement for enough thermal ionization to initiate a robust MRI is a 
central temperature T > 1000 K (Gammie 1996), comparable to our observationally-based 
estimates (Te//~800 K at i? ~ 200 Rq). 
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5.3. Thermal and other instabihties 



The thermal instability model was originally suggested to account for outbursts in dwarf 
nova systems (Faul kner, Lin, &: Pap aloizou 1983) and has attractive properties for explaining 
FU Ori outbursts (IBell &: Linlll994l ). However, the observed large outer radius, 0.5 — 1 AU, 
of the high accretion rate portion of the FU Ori disk poses difficulties f or th e pure thermal 
instability model. Even in the situation explored by iLodato fc Clarkd (120041 ). in which the 
thermal instability is triggered by a massive pla net in the inner disk, the outer radius of the 
high state is still the same as in lBell fc LinI (119941 ). ~ 20 R©. The relatively high temperatures 
(2 — 4 X 10^ K) required to trigger the thermal instability (due to the ionization of hydrogen) 
are difficult to achieve at large radii; they require large surface densities which in turn produce 
large optical depths, trapping the radiation internally and making the central temperature 
much higher than the surface (effective) temperature. 

To illustrate the problem, we calculated thermal equilibrium "S" curves (e.g., Faulkner 
et al. 1983) for our disk parameters at R = 210Rq but various values of a. The equilibrium 
curve represents energy balance between viscous energy generation and radiative losses, 



rad 



(6) 



Thus, 



4 



(7) 



where S is the surface density and Cg is the sound speed at the local central disk. For 
an optic ally-thick disk , we ha'v 
opacity ( Hubenylll99ol ). Then, 



an optic ally-thick disk , we have T^=|r/jTgjj, where tr=krT, and k is the Rosseland mean 



27 



where TZc is the gas constant, and fi is the mean molecular weight . As long as we know 
Tc, we can derive Tg// and the corresponding M at this radius. Because the Rosseland 
mean opacity is also dependent on temperature and pressure, iteration of the disk structure 
calculation is needed to derive the equilibrium curve. 

We show in Figure fTOl the equilibrium curve at R=210 R© for five values of a: 10~^, 
10^^, 10"^, 10"'^, 10"^. In discussing our results, we define to be the highest stable low 
state for a given a. If at any radius the surface density increases above S^, the disk can 
no longer stay stable on the lower branch and strong local heating begins. If the a of the 
low state is the same as the a of the high state, then a>10~^; in this case, to trigger the 
thermal instability at S^, the mass accretion rate in the low state would need to be ~10~^ 
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(see Figure [To]), which is higher than the M of the high state that we have determined. To 
have a mass accretion in the low state at least one order of magnitude lower than that of the 
high state, we would n eed a<10~ -^ for th e low state (see Figure [TU]). This is comparable to 



the low state a used by lBell &: Lid (119941 ). but very much lower than our high state estimate. 
Changing the outer radius to R~100 R0 does not change these results by more than a factor 
of 10. 

With a=10~'^, the critical surface density at 210 is ~ 10^ g cm~^. Then we can 
estimate the Toomre Q parameter, 

Q = ^-o.i, (9) 

where Q < 1 implies gravitational instability. Thus, a thermal instability model for FU Ori 
implies a massive, probably gravitationally-unstable disk, which raises the question: could 
the outburst be driven by gravity instead? 



Armitage et al.l (120011 ) proposed an alternative model of disk outbursts i n which rnateria l 



piles up at radii of order 1-2 AU, consistent with layered disk models ( iGammid Il996l ). 
achieving repeated outbursts with high accretion on size scales much closer to our inferred 
Rout- In this model, gravitational torques lead to increased accretion, which in turn heats up 
the inner disk until temperatures are high enough (800 K in their mo del) to trigger th e MRI, 



which results in very much higher viscosities and rapid accretion. The lArmitage et al.l (120011 ) 



model does not give exactly FU Ori-type outbursts; instead, the accretio n events last for 



10^ yr and achieve peak accretion rates of only ~ 10 ^Mq yr ^. However, lArmitage et al. 



( I2OOII ) suggested that thermal instabilities might be triggered by their gravitationally-driven 
outbursts. Simple models with somewhat different parameters demonstrate that it is indeed 
possible to trigger thermal instabilities in the context of this general model (Gammie, Book, 
& Hartmann 2007, in preparation), resulting in outbursts more similar to that of FU Ori. 

Vorobyov & Basu (2005,2006) suggest FU Ori outburst could be explained by the accre- 
tion of clumps formed in a gravitationally unstable disk. The disk becomes gravitationally 
unstable because of the infall of the matter from the envelope. The spiral arms and clumps 
form and grow in the unstable disk. The redistribution of mass and angular momentum 
leads to centrifugal disbalance in the disk and further triggers the outburst when the dense 
clumps are driven into the central star. The authors derive a mass accretion rate at 10 AU 
of ~10~^MQ/yr, which could result in too much emission from the outer disk; this issue 
requires further exploration. 

The rapid rise time of the outburst of FU Ori in the B band, ~ 1 year, provides the 
best evidence for a thermal instability. We can test this idea very crudely by using our 
steady disk model to estimate the range of radii in the disk which can contribute to the B 
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band magnitude in the light curve, finding that only the inner R < 2OR0 disk are involved 
If we then assumed that the rise in B light is due to the inward propagation of a thermal 



insta bility front across this region at a (maximum) speed ~ aCg (ILin et al.lll985l : iBell et al. 



19951 ). and assume a sound speed characteristic of 2 — 4 x 10^ K (the minimum temperature 
for instability), then we derive a ~ 0.1. While this calculation is very crude, it does suggest 
consistency with the idea of a thermal instability with a viscosity parameter comparable to 
that which we inferred from the decaying light curve (§5.2). 

Finally, we estimate the amount of mass in the inner disk. An upper mass limit for 
the inner disk can be estimated by assuming Q ~ 1 (the approximate limit for strong 
gravitational instability) at i? ~ 200 R0. Adopting a central temperature T ~ 1000 K, the 
gravitational instability limit leads to an approximate surface density limit S < 10^gcm~^ 
and thus a maximum inner disk mass ~ vrSi?^ ~ 0.05Mq. This leads us to an inner disk mass 
estimate O.O5M0 > Mjn > O.OIM©, where the lower limit comes from requiring the outburst 
mass accretion rate of ~ lO~^M0yr~^ to be sustained for the observed decay timescale of 
~ 100 yr. 



6. Conclusions 

Using the latest opacities and the ODF method, we have constructed a new detailed 
radiative transfer disk model which reproduces the main features of the FU Ori spectrum 
from ~4000yl to about 8 fim. The SED of FU Ori can be fitted with a steady disk model and 
no boundary layer emission if Ay~1.5, a somewhat lower valued extinction than previously 
estimated. A larger Ay would imply extra heating at small radii which would result from 
dissipation of kinetic energy as material accretes onto the central star. With Av=1.5, the 
inclination estimated from near-IR interferometry and the observed rotation of the inner 
disk, we estimate that the central star has a mass ~ O.SMq (typical of low-mass TTS) and 
a disk mass accretion rate of ~ 2.4 x 10~^Mq yr^-*^. 

Relying on the Infrared Spectrogr aph Spectrum of FU Ori from the Spitzer Space Tele- 



scope presented by lGreen et al.l (120061 ) . we estimate that the outer radius of the hot, rapidly 
accreting region of the inner disk is ~1 AU, although the irradiation effect may decrease 
this value by no more than a factor of two(Zhu et al., in preparatio n) . Either way, this 



is inconsistent with the pure thermal instability models of iBell fc LinI (jl994j ) who adopted 
a low a viscosity parameter ~10~'^-10~^. If we assume that the observed decay timescale 
of FU Ori (~100 years) is the viscous transport time from the outer edge of the hot re- 
gion, then we derive a ~ 0.2 — 0.02, , comparable to the predictions of simulations of the 
magnetohydrodynamic (MHD) turbulence. The effective temperature at the outer edge of 
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the hot region in our model is ~800 K ; the presence of absorption features shows that the 
central disk temperature must be higher. Thus, our model suggests that the central disk 
temperature is greater than the Tc~1000 K required to maintain full MHD turbulence by 
thermal ionization out to ~ 1 AU. 

We show that pure thermal instability models have difficulty in explaining the outburst 
of FU Ori, and suggest that models including gravitational torques and MRI activation, such 



as those of lArmitage et al.l (120011 ). are more promising. 
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Table 1: Best fit model 
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Table 2: Optical lines used to determine the spectral type 



Lines ^ 


Metals 


Spectral type 


4047.00 


Fe I+Sc I 


F8 


4226.00 


Ca 1 


Gl 


4271.00 


Fe 1 


G4 


4305.00 


CH(Gband) 


Gl 


4458.00 


Mn I+Fe I 
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5329.00 


Fe I 


G2 


5404.00 


Fe I 


G3 


6162.00 


Ca I+TiO 


F9 



"From (Hernandez 



( 20051) 
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Fig. 1. — The spectral energy distribution of FU Ori. Collected photometric and spectro- 
scopic data of FU Ori with A^=1.5 (see §4). The solid circles with error bars are the UBVR 
photometry from the Maidanak Observatory, the squares are 2MASS JHK photometry, and 
the open circles are JHK'L' photometry from lReipurth &: AspinI (120041 ) . The dashed curve is 
the FAST spectrum from 3650 A to 7500 A , the dotted curve is the KSPEC spectrum from 
1.15 yum to 2.42 u m (IGreene fc Ladal 119961) . and the solid curve is the SpeX spectrum from 
2.1 fim to 4.8 fim (iMuzeroUe et al.ll2003l ). both scaled to the 2M ASS photome t ry. F inally, 
the dot-dash curve is the IRS spectrum from 5 fim to 35 fim from lGreen et al.l (120061 ). The 
spectra were put on as absolute flux scale as discussed in §2. 
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Fig. 2. — Our calculated Rosseland mean opacity as a function of temperature for so- 
lar composition. Each curve is labeled with the value of log(p/T|) to be compared with 
Alexander fc Ferguson! (119941 ). where Tg is the temperature in millions of degrees. 
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Fig. 3. — The spectrum predicted by the a nnulus with T^f f= 5300 K (sohd hne) and the 
spectrum of SAO 21446 which is a Gl I star (jjacoby et al.lll984j ) ( dotted hne). The spectra 
are scaled to the same flux at 7000 A. 
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Fig. 4. — Emergent intensities from selected annuli of the disk model calculated for a 55 
degree emergent angle. Effective temperatures of the annuli at 1, 3, 5, 6, 8, 10, 15, 20, 
25, and 30 Rj are 6420, 4660, 3400, 2480, 2130, 1630, 1310, 1110, and 977 K, respectively. 
Molecular absorption features in the near-IR arise from cool regions of the disk, which do 
not contribute to the optical spectrum. 
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Fig. 5. — Observed UBVR photometry and optical spectra with different Ay: 2.2,1.9,1.5,1.3 
(dotted Unes, from top to bottom) and model spectra with Tj^ax' 7240, 6770, 6420, 5840 
(solid lines, from top to bottom). The total flux spectra are produced by weighting the 
intensities (Figure 4) by the appropriate annular surface areas. 
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Fig. 6. — Comparison between observed data and disk models for different values of the 
outer radii. All the black lines and dots are observed spectra and photometry data as shown 
in figure [H All the color lines are model spectra. The parameters of the disk models are 
inner disk radius 5 R©, and MM~ 7.2xlO^^MQ/yr. The corresponding Tmax is 6420 K. 
The four model spectra from top to bottom are for the disk models truncated at 5000 R©, 
210 Rq, 50 Rq and 20 Rq. The solid lines are model spectra for a high mass accretion rate 
region with outer radii as shown. The dotted lines are spectra for models where an outer 
region with mass accretion one order of magnitude smaller has been added. For the disk 
model truncated at 5000 Rq, the dotted line is coincident with the solid line. The extinction 
parameter is Ay =1.5. 
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Fig. 7. — Two truncated models fitting the observations corrected by Av=2.2 and Ay=1.5. 
Only optical photometries and IRS spectra (5-30 microns) are shown. The square points 
are the photometries corrected by Ay =2. 2, while the round points are the photometries 
corrected by Ay =1.5. The IRS spectra corrected by these two Ay are coincident. The solid 
line is the model fitting to observations with Ay=2.2, while the dotted hne is the model 
fitting to observations with Ay =1.5. Both of the two models have the same outer radius of 
the high state Ro„t=210 Rq. 
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Fig. 8. — CO first overtone absorption ban d around 2.3 fim. The s olid line is the observed 
CO first overtone absorption feature from iHartmann et al.l (120041 ) and the dotted line is 
the model spectrum after being convolved to the same spectrum resolution. The turbulent 
velocity is 4 kms~^ and the inclination angle is 55 degrees. 
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Fig. 9. — The mid-infrared part of the observed spectra in Figure [1] compared with the model 
spectra. The sohd hne is the model with Ro„t=210 R0 and the dotted line is the model with 
1 % dust opacity (dust opacity reduced by a factor of 100) and RoMt=225 R©. We do not 
have the CO fundamental opacities in our model and thus cannot reproduce these strong 
spectral features at 4.6 /xm. 
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Fig. 10. — The equilibrium curves for thermal-instability at R=210 Rq with five values of a: 
10-^10-^10-^10-^10-^(from top to bottom). The dotted line is M = 2.4 x lO-^M^/yr, 
set to be one order of magnitude lower accretion rate than the accretion rate of the inner 
disk. If the outer disk M were one order of magnitude higher than this (so that the outer 
disk and inner disk accretion rates are the same) the mid-infrared flux would be much larger 
than observed (see discussion in §4). The open circles denote the critical surface densities 
T,A for triggering the thermal instability at each value of a (see §5.3) 



